Charmonium in strongly coupled Quark-Gluon Plasma 
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The growing consensus that a strongly-coupled quark-gluon plasma (sQGP) has been observed at 
the SPS and RHIC experiments suggests a different framework for examining heavy quark dynamics. 
We present both semi-analytical treatment of Fokker-Planck (FP) evolution in pedagogical examples 
and numerical Langevin simulations of evolving cc-pairs on top of a hydrodynamically expanding 
fireball. In this way, we may conclude that the survival probability of bound charmonia states is 
greater than previously estimated, as the spatial equilibration of pairs proceeds through a "slowly 
dissolving lump" stage related to the pair interaction. 

PACS numbers: 



I. INTRODUCTION 



Overview 



Charmonium suppression is one of the classical probes 
used in heavy ion collisions. Since charm quark pairs 
originate during early hard processes, they go through all 
stages of the evolution of the system. A small fraction 
of such pairs ^ 0(10^^) produces bound cc states. By 
comparing the yield of these states in AA collisions to 
that in pp collisions (where matter is absent) one can 
observe their survival probability, giving us important 
information about the properties of the medium. 

Many mechanisms of J/^p suppression in matter were 
proposed over the years. The first, suggested by one of 
us in 1978 [l|, is (i) a gluonic analog to "photo-effect" 
gJ/i^ cc. Perturbative calculations of its rate Q pre- 
dict a rather large excitation rate of the charmonium 
ground state. Since charmonia are surrounded by many 
gluons in QGP, this leads to a conclusion by Kharzeev 
and Satzfs] that nearly all charmonium states at RHIC 
should be rapidly destroyed. If so, the observed J/ip 
would come mostly from recombined charm quarks at 
chemical freezout, as advocated in (3]. 

However this argument is only valid in weakly-coupled 
QGP, in which the charm quarks would fly away from 
each other as soon as enough energy is available. As 
we will show below, in strongly-conpled QGP (sQGP) 
propagation of charmed quark is in fact very different. 
Multiple momentum exchanges with matter will lead to 
rapid equilibration in momentum space, while motion in 
position space is slow and diffusive in nature. Persis- 
tent attraction between quarks makes the possibility of 
returning back to the ground state for the J/ ip quite sub- 
stantial, leading to a substantial survival probablity even 
after several fm/c time in sQGP. 

Another idea (ii) proposed by Matsui and Satz [5] fo- 
cuses on the question of whether charmonium states do 
or do not survive as a bound state. They argued that be- 
cause of the deconfinement and the Debye screening, the 
effective cc attraction in QGP is simply too small to hold 



them together as bound states. Quantum-mechanical 
calculations by Karsch et al Q and others have used the 
free energy, obtained from the lattice, as an effective po- 
tential (at T > Ta) 



F{T, r) 



3r 



exp(-MD(r)r) + F{T, oo) (1) 



They have argued that as the Debye screening ra- 
dius decreases with T and becomes smaller than 
the root mean square radii of corresponding states 
X, •^/V'l "^"i Tf'j T..., those states should subsequently 
melt. Furthermore, it was found that for J/^p the melt- 
ing point is nearly exactly Tc, making it a famous "QGP 
signal" . 

These arguments are correct asymptotically at high 
enough T, but the central issue is what happens at (so 
far experimentally accessible at RHIC) T = (1 — 2)Tc. 
Dedicated lattice studies extracted quarkonia spectral 
densities using the so called maximal entropy method 
(MEM) to analyze the temporal (Euclidean time) cor- 
relators. Contrary to the above-mentioned predictions, 
the peaks corresponding to rjcJ/ili states remains ba- 
sically unchanged with T in this region, indicating the 
dissolution temperature is as high as « (2.5 — ?))T(.. 
Mocsy et al J,3] have used the Schrodinger equation for 
the Green function in order to find an effective potential 
which would describe best not only the lowest s-wave 
states, but the whole spectral density. Recently they 
have argued that a near-threshold enhancement is hard 
to distinguish from a true bound state: according to these 
authors, the above mentioned MEM dissolution temper- 
ature is perhaps too high. 

Another approach to charmonium in heavy-ion colli- 
sions, taken by Grandchamp and Rapp Q, does not rely 
on the perturbative calculation of the excitation cross- 
section. Charm rescattering is enhanced by formation 
of bound states in QGP. J/ij] lifetimes were calculated 
at various temperatures using heavy quark effective the- 
ory [7]: the resulting widths are typically a few hundred 
MeV. If so, the total number of rescatterings of J/ip in 
the fireball during its lifetime 10 fm/c) is large (10-30). 
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This model still has fairly large cross-sections for J /ip - 
annihilation, so in their so-called two-component model, 
many of the final charmonia measured are required to 
originate from statistical coalescence of single charm in 
the plasma into charmonium states. 

There are also other quantum-mechanical studies of 
the issue since the pioneering paper Q. Zahed and 
Shuryak [11] argued that one should not use the free 
energy F[T, r) as the effective potential, because it cor- 
responds to a static situation in which infinite time is 
available for a "heat exchange" with the medium. In 
the dynamical real-time situation they proposed to think 
in terms of level crossing and Landau-Zener formalism, 
widely used in various quantum-mechanical applications. 
In the "fast" limit (opposite to the "adiabatically slow" 
one) all level crossings should instead be ignored. This 
corresponds to retaining pure states (described by a wave 
function rather than density matrix) without the entropy 
term in F, which is nothing else but the internal energy 

V{T, r) = F(T, r) - TdF{T, r) /dT ^ F + TS (2) 

instead, as an alternative effective potential. Such poten- 
tial V{T, r) (extracted from the same lattice data) leads 
to much more stable bound states, putting charmonium 
melting temperature to higher T ~ 3Tc. A number of au- 
thors [l2| have used effective potentials in between those 
two limiting cases. However, as it will be clear from what 
follows, we think it is not the bound states themselves 
which are important, but kinetics of transitions between 
them. In a nutshell, the main issue is how small is the 
separation in the cc pair when the QGP is over, not in 
which particular states they have been during this time. 

The heavy QQ potential depends not only on the tem- 
perature but also on the velocity of the cc pair relative to 
the medium. This effect has been studied e.g. by means 
of the AdS/CFT correspondence in [l^ and it was found 
that the bound state should not exist above a certain 
critical velocity. So, if the existence of a bound state is 
truly a prerequisite for J/tp survival, one would expect 
additional suppression at large pt- This goes contrary 
to the well-known formation-time argument fl&l and the 
experimental evidence, indicating the disappearance of 
the J/ip suppression at large pt ■ We think this as a good 
example of how important the real-time dynamics of the 
cc pair in medium is: and indeed we are going to follow 
it below from its birth to its ultimate fate. 

Let us now briefly review the experimental situation. 
For a long time it was dominated by the SPS experiments 
NA38/50/60, who have observed both "normal" nuclear 
absorption and an "anomalous" suppression, maximal in 
central Pb+Pb collisions Since at RHIC QGP has 
a longer lifetime and reaches a higher energy density, 
straightforward extrapolations of the naive J/ip melting 
scenarios predicted near-total suppression. And yet, the 
first RHIC data apparently indicate a survival probabil- 
ity similar to that observed at the SPS. 

One possible explanation ^20] is that the J/ip suppres- 
sion is cancelled by a recombination process from unre- 



lated (or non-diagonal) cc pairs floating in the medium. 
However this scenario needs quite accurate fine-tuning 
of two mechanisms. It also would require rapidity and 
momentum distributions of the J/ip at RHIC to be com- 
pletely different from those in a single hard production. 

Another logical possibility [2l[ is that the J/ip actually 
does survive both at SPS and RHIC: while the anomalous 
suppression observed may simply be due to suppression of 
higher charmonium states, ip and x, feeding down about 
40% of J/ip in pp collisions. These authors however have 
not attempted to explain why J/ip survival probablity 
can be close to one. 

This is precisely the goal of the present work, in which 
we study dynamically how survival of J/ip happens. We 
will see that it is enhanced by two famous signatures of 
sQGP, namely (i) a very small charm diffusion constant 
and (ii) strong mutual attraction between charmed quarks 
in the QGP. We found that J/ip survival through the 
duration of the QGP phase r ~ 5fm/c is about a half. 

The sequence of events can be schematically described 
as a four-stage process 



(cc production) 



/i. 



itial 



(mom. relaxation) 



^quasi— cquilibriuin(3) 



(leakage^ 



/f 



final 



(projection) 



j/ib 



A new element here is a two-time-scale evolution, includ- 
ing rather rapid momentum relaxation to a quasiequilib- 
rium distribution which differs from the equilibrium one 
at large enough distances. 



B. Charmonium potentials at T > Tc 

The interaction of the cc pair will play a significant 
role in this paper, and thus we briefly review what is 
known about them. The details can be found in the orig- 
inal lattice results: we will point out only the most im- 
portant qualitative features. 

Perturbatively, at high T one expects a Coulomb-like 
force, attractive in the color singlet and repulsive in the 
color octet channel, with the relative strengths 8:(-l) (so 
that color average will produce zero effect). As shown 
by one of us many years ago[lj, the Coulomb forces are 
screened by the gluon polarization operator at distances 

-i/gT.^ 

Quantitative knowledge of the interaction comes from 
large set of lattice measurements of the free energies asso- 
ciated with a pair of heavy quarks in an equilibrium heat 
bath. These data include both results by the Bielefeld- 
BNL group and in dynamical QCD with Nf 2 by Aarts 
et al. [3. 

At T > Tc one is in a deconfined phase, so at large 
quark separations one expects effective potentials to go 
to a finite VooiT). Yet when the value of this potential 
significantly exceeds the temperature, the actual proba- 
blity of quark separation is small ~ exp(— Voo(r)/T). As 
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FIG. 1: (Color online.) Parametrization of the potential for 
lattice data from [23|. From top to bottom: the potential at 
T/Tc = 1.02, 1.07, 1.13, 1.18, 1.64 . 



we already mentioned in the preceding section, the appro- 
priate potential for dynamical cc pair is not yet definitely 
determined, with suggestions ranging from free energy to 
potential energy measured on the lattice. 

The difference between the two - the entropy associ- 
ated with cc pair - is very large near Tc, reaching the 
value S* ~ 20 at its peak, ft means that a huge number 
of excited states ~ exp(S') ~ exp(20) would be excited 
in adiabatically slow motion of the pair. We think that 
in realistic motion of cc pair much less states are ac- 
tually excited, and thus, following we will use the 
potential energy instead of free energy. A cost of this is 
much larger potential barrier, reducing cc dissociation. 
Indeed, Voo{T) near Tc is large, reaching about 4 GeV at 
its peak. 

For the simulation we need a parametrization of the 
heavy quark- antiquark potential above deconfinement as 
a function of temperature and separation. We use the 
same parametrization as in [is!]: 



y(r,T) 



(4) 



and fit to quenched-QCD lattice data in a temperature 
range of 1.02Tc to 1.64Tc, assuming a is constant in tem- 
perature, /i varies linearly with temperature over this 
range, and C{T) cx {T/Tc - 0.98)~\ so that it peaks 
sharply at Tc as C/jnf does. The result shown in Fig llBT a) 
is for 



m(T) = (0.03 + 0.006T/Tc) GeV^ 



cr = 0.22 GeV^ 



C(T) 



O.lSGeV 
T/Tc ~ 0.98 



(5) 



(6) 



(7) 



with a set to j^. 

As you can see from Fig ll Bf a). this fit is not perfect. 
While it is easy to fit a function to a single temperature's 
data set, it is hard to find an adequate fit across temper- 
atures. This fit however will prove sufhcient, especially 
since it is relatively good for the small separations of in- 
terest. 

The classical Boltzmann factor exp(J7(r)/T) for a 
Coulomb potential leads to a non-normalizable distribu- 
tion. Quantum mechanics prevents this, which can be 
crudely modeled by an effective potential 



Veff = h^/Mr^ + V{T,r) 



(8) 



which includes the so called localization energy. Dusling 
and one of us have determined more accurate effec- 
tive quantum potentials, following ideas of Kelbg and 
others and performing path integrals. Perhaps those 
should be used in future more sophisticated simulations. 
In the present simulations we simply turn off the force 
below 0.2 fm, approximating the effective potential by a 
constant. 

During the simulation, we allow our pairs to exist as 
either color singlets or color octets. While in a zero- 
temperature pp collision, confinement completely sup- 
presses a cc pair's probability of existing in a color octet 
state, in the deconfined phase this possibility must be 
considered. We initially create a 1:8 ratio of color singlet 
to color octet pairs, as is expected statistically, and then 
during each timestep we decide if the pair exists as a 
color octet by comparing a random number between zero 
and one to ^ exp(— ([/g — Ui)/T), with the color singlet 
and octet potential energies determined from [3^. In 
other words, the pairs exist in thermal equilibrium in 
color space at all times. When a pair exists as a color 
octet, they do not interact in our simulation as the spatial 
variation of the octet free energies is quite small. 

Naively one would expect this to create a large dif- 
ference with results where color singlet is assumed for 
all pairs, however this is not the case. One sees this by 
noticing that for the temperatures and distances of our 
interest (distances of separation for pairs likely to go into 
the J/'0-state), ^^^^^ ^ 10, meaning that the octet state 
is suppressed by orders of magnitude. 

Lattice potentials do not follow the simple perturbative 
relation Vi = — SVg between singlet and octet potentials 
mentioned above. While the color singlet channel dis- 
plays a significant attraction, the color octet channel has 
a potential which is remarkably flat (r-independent). On 
the basis of this feature of the lattice data we will ignore 
the force in the octet channel in the simulation. 



C. Charm diffusion constant 

Loosely speaking, the effect we are after in this work 
can be express as follows: the medium is trying to pre- 
vent the outgoing quarks from being separated. It has 
been conjectured by one of us that charm quark pairs 
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should get stopped in QGP ^25*1. RHIC single electron 
data suggest that charm quarks, and probably b quarks 
as well, are indeed equilibrating much more effectively 
than it was thought before. 

The first question one should address is how a charm- 
anticharm pair moves in sQGP, and what is the proba- 
blity to find the cc in close proximity after time t. Charm 
quarks are subject of three forces: (i) the drag force, try- 
ing to reduce the difference between quark momentum 
distribution and that in (local) equilibrium quark and 
matter velocities 

(ii) stochastic (or Langevin) force from a heat bath, lead- 
ing to thermal equilibration in momenta and spatial dif- 
fusion of the charm quarks 

(iii) the cc mutual interaction. For pedagogical reasons it 
is useful to include them subsequently, first for a station- 
ary non-floating matter at fixed temperature (using the 
Fokker-Planck equation) and then for a realistic nuclear 
geometry with a hydrodynamically expanding fireball. 

But before we do so, let us briefly remind why the case 
of heavy (charmed) quarks is so special. (More detailed 
discussion of that can be found in Moore and Teaney 
[1^.) A collision of a heavy quark with quasiparticles of 
QGP leads to change in its momentum Apng ~ T, so 
that the velocity is changed by a small amount 

AvHQ - T/Mhq « 1 (9) 

Therefore the velocity of the charm quark can only 
change significantly as a result of multiple collisions, in 
small steps. Thus the process can be described via ap- 
propriate differential equations such as the Fokker-Planck 
or the Langevin equations. Similarly one can argue that 
spatial diffusion of a heavy particle can also be described 
in this way, because the change in position between col- 
lisions are small and uncorrelated. 

An assumption necessary for Langevin dynamics to 
hold is that the "kicks" are random (uncorrelated). As 
explained above, for a single heavy quark it follows from 
the inequality M >> T, which guarantees that the quark 
relaxation time is long compared to the correlation time 
in matter. For a cc pair we need an additional require- 
ment, that random forces on each quark can be treated 
as mutually uncorrelated. In order to see how good this 
approximation is, one should compare the spatial (equal 
time) correlation length ^(T) in the medium to the typ- 
ical distance between quarks for paths which eventually 
(at the end of the plasma era) will become charmonia. 
We will provide two different estimates for the former 
quantity, which view QGP either as a perturbative "gas" 
or a strongly coupled "liquid" . 

In a perturbative gas of gluons, the mass is small and 
in the lowest order momentum distribution is thermal. 
Thus the maximum of the momentum distribution is at 
p w 2.7T, about 1 GeV at the initial RHIC condition. 
The corresponding half wavelength (the region where 
the field keeps at least its sign) 



In the liquid regime quasiparticles do not successfully 
model the degrees of freedom. However we do have 
phenomenological information about spatial correlations 
from hydrodynamics, which propose the so called "sound 
absorption length" as a measure above which different 
matter "cells" decorrelate. It is 



with rj/s, the dimensionless ratio of viscosity to entropy 
density. Empirically, RHIC data are well described if it 
is of the order of I/Stt (the AdS/CFT strong-coupling 
limit), which suggests a spatial correlation length one 
order of magnitude smaller ^(T = .4 GeV) ^ 0.05 fm. 

Since the distance between c and c which eventually 
become J/^ is about .5 fm, and since there are good rea- 
sons to believe the latter "liquid" estimate of ^ is closer 
to reality, we conclude that the main Langevin assump- 
tion - the independence of random forces for c and c - is 
well justified. Furthermore, one may think that the same 
assumption would even hold for b and &, although with 
worse accuracy. 

Since this small parameter in Eq[S] is central to what 
follows, let us remind for orientation of the reader that 
at RHIC we speak about the ratio for charmed quark 
T/Mc = (1/6 - 1/5), or T/Mt = (1/20 - 1/15) for b 
quark. 

Although RHIC experiments with charm quarks in- 
clude direct reconstruction of charmed mesons D, D* by 
the STAR collaboration, so far the existing vertex detec- 
tors are not sufficient for doing it effectively (upgrades 
are coming). Therefore the most relevant data on charm 
is based on observation of single electrons from heavy 
quark weak semileptonic decays. Apart of electromag- 
netic backgrounds, we do not really know whether elec- 
trons come from c or 6 decays: it is believed (but not 
yet proven) that the boundary between two regimes is at 
Pf w 4 GeV. Two experimentally observable quantities 
are (i) the charm suppression relative to the parton model 
(no matter) R\j^ and (ii) the azimuthal asymmetry of the 
electrons relative to impact parameter t;| =< cos{2(j)) >. 

Several theoretical groups have analyzed these data, in 
particular Moore and Teaney [26j provided information 
about the diffusion constant of a charm quark Dc by 
Langevin simulations. A conclusion following from this 
work is that both i?^^ and wf observed at RHIC can be 
described by one value for the charm diffusion constant 
in the range 

Z^c (27rr) = 1.5 - 3. (12) 

This can be compared with the perturbative (collisional) 
result at small Us 

jjpQCD (27rT) = 1.5/al. (13) 
Assuming that the perturbative domain starts^^ some- 



A/2 = 7r/pw .6fm- ^(r= .4GeV) (10) 



Recall that at A/3as = 1/2 two scalar quarks should fall towards 



5 



where at tts < 1/3 one concludes that the empirical value 
is an order of magnitude smaller than the perturba- 
tive value. 

There are two studies of the diffusion constant at 
strong coupling. One comes from AdS/CFT correspon- 
dence Q (and thus the results are for A/'=4 supersym- 
metric Yang-Mills theory): the final expression found by 
Casalderrey-Solana and Teaney 28] is 

Dhq = % (14) 

It is nicely consistent (via the Einstein relation) with the 
calculated drag force 

dP 



(15) 



by Herzog et al [29|. One assumption of this calculation 
is that the 't Hooft coupling is large g'^Nc ^ 1, which 
means the diffusion constant is parametrically small, 
much less than momentum diffusion in the same theory 
Dp = r]/{e+ p) ^ 1/AttT. This result is also valid only 
for quarks heavy enough 



M 



(16) 



which only marginally holds for charm quarks. 

Let us see what these numbers mean for RHIC (as- 
suming that they are valid for QCD). The 'tHooft cou- 
pling g'^Nc — as^TrNc ~ 20 — 40 is indeed large, while 
Meff ~ 1 — 2 GeV is not really small as compared to the 
charm quark mass: thus the derivation is only marginally 
true. Yet we proceed and get 



{2ttT) = 4/vV^ - 0.5 - 1, 



(17) 



which is in right ballpark of phcnomcnological numbers. 

Another approach to transport properties of a strongly 
coupled plasmas is classical molecular dynamics. A clas- 
sical non-Abelian model for sQGP was suggested by Gel- 
man et al ^3^, and recently Liao and Shuryak "sij have 
added magnetically charged quasiparticles. Those cal- 
culations also find that the diffusion constant D rapidly 
decreases as a function of the dimensionless coupling con- 
stant r as a power 



D . 



0.6-0.8 



(18) 



in a liquid domain F = 1 — 100. Qualitatively it is similar 
to the \/y/gWc of the AdS/CFT result Unfortu- 
nately, at this moment there is no deep understanding of 
the underlying mechanism of both strong coupling calcu- 
lations. 

But this is well beyond the aims of the present work: 
in what follows we will use Dc {2nT) = 1.5, 3 as a range 
for our best current guess. 



each other, according to the Klein-Gordon equation: so this is 
clearly not a perturbative region. 



II. CHARMONIUM IN SQGP: 
FOKKER-PLANCK FORMALISM 

Before we described realistic Langevin simulations of 
RHIC collisions it is useful to see first the basic features 
analytically. In this section we describe how (1) the high 
drag coefficient rjc causes rapid thermalization of the ini- 
tially hard momentum distribution from PYTHIA, and 
(2) the small diffusion coefficient Dc, inversely propor- 
tional to r/c, combined with (3) the use of the static QQ 
internal energy instead of the free energy for the pair in- 
teraction leads to a slow dissolution of a large peak in 
the distribution in position space, causing less suppres- 
sion than expected by more naive models. 

To see the first feature, rapid thermalization in mo- 
mentum space, consider a cc pair with the quarks emerg- 
ing back-to-back from a hard process with momentum 
p >> T, therefore having relative momentum Ap 2p. 
For p >> T, as is the case initially for charm quarks 
created at the RHIC, we may neglect the random kicks 
that the quarks experience from the medium and consider 
only the drag: = -VcP, where r]c = j^^- This 

leads to a very simple formula for the relative momentum 
of a pair at early times: 



Ap(t) = exp(-?7,T)Ap(0) 



(19) 



Using the AdS/CFT value for the diffusion coefficient, 
this leads to a drag coefficient rjc ~ .6 GeV. 

Let us now examine the distribution of cc pairs cre- 
ated with PYTHIA pp-event generation at energies of 
200 GeV (see more on this in Section [Ill| . The initial 
transverse momentum distribution ^ is broad compared 
to the thermal distribution, and therefore we may ap- 
ply Eq. [Tni for early times. We parametrize the initial 
relative momentum distribution as in Section HITB] and 
replace Ap with Apexp(77cT), which is the formula for 
a pair's initial relative momentum in terms of its rel- 
ative momentum at proper time t. Next, we compute 
the overlap of this distribution at early times with the 
J/ip -state's Wigner quasi-probability distribution to de- 
termine the probability of a random cc taken from this 
ensemble to be measured as a J/ip particle, an approach 
detailed in Appendix [Bl 

Fig. |TT] shows this probability as a function of time 
for rjD — 0.88. The initial value of the projection is 
0.8%, on the same order of magnitude as the experimen- 
tal value of 1% for ^ obtained from [11], [33]. The 
projection increases as the PYTHIA-with-drag distribu- 
tion narrows but then drops when the width shrinks more 
than the J/ip^s width in momentum space and as the 



^ The initial rapidity distribution is wide also, but since longitu- 
dinally Bjorken-like hydrodynamics starts immediately, a charm 
quark finds itself with comoving matter with the same rapidity, 
and thus has little longitudinal drag. Transverse flow is slow to 
be developed, thus there is transverse drag. 
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FIG. 2: Probability of a cc pair going into the J/tp-state vs. 
time, for very early time. 



quark pair's separation increases. Once the time reaches 
about 1 fm/c, the probabihty for a pair to go into a 
J/ip-state is about where we started, and we quit look- 
ing at this approach after this much time because the 
mean transverse momentum for a quark is the thermal 
average. 

After this first 1 fm/c of the QGP phase, the 
cc distribution has thernialized in momentum space and 
the evolution in position space (diffusion) needs to be 
examined. The root mean square distance for diffusive 
motion is given by the standard expression 



(20) 



where r is the proper time and the interaction between 
the quarks has been neglected. The "correlation volume" 
in which one finds a quark after time r is 

4-7r 

Korr = —{GD.rf/' (21) 

and one may estimate for the probability of the cc pair 
to be measured in the J/4' -state as 



(22) 



So neglecting the pair's interaction leads to a small prob- 
ability that J/^p -states will survive by the hadronization 
time at the RHIC (r ~ 10 fm/c), even for small values of 
the diffusion coefficient. 

To get an idea for how this simple result is changed by 
the inclusion of an interaction between the constituent 
quarks in a given cc -pair, let us examine the Fokker- 
Planck equation for the cc distribution in relative posi- 
tion: 

dP 

'dt 



or or 



FIG. 3: (Color online.) Numerical solution of the 
one-dimensional Fokker-Planck equation for an interacting 
cc pair. The relaxation of the initial narrow Gaussian distri- 
bution is shown by curves (black, red, brown, green, blue, or top 
to bottom at r=0) corresponding to times t = 0, 1,5, 10 fm, 
respectively. 



(23) 



where /o(f) oc exp{—Vef[ (r) /T) is the equilibrium dis- 
tribution in the magnitude of relative position r. By 
substituting the potential shown above at T = 1.25rc 
and Dc x {2ttT) = 1 into the Fokker-Planck equation 
(for demonstration in a single spatial dimension only) we 
solve it numerically and find how the relaxation process 
proceeds. A sample of such calculations is shown in Fig. 
im It displays two important features of the relaxation 
process: 

(i) during a quite short time T, 1 fm the initial distribu- 
tion (peaked at zero distance) relaxes locally to the near- 
equilibrium distribution with two peaks, corresponding 
to optimal distances of the equilibrium distribution 
where the effective potential is most attractive; 

(ii) the second stage displays a slow "leakage", during 
which the maximum is decreasing while the tail of the 
distribution at large distances grows. It is slow because 
the right-hand side of the Fokker-Planck equation is close 
to zero, as the distribution is nearly /q. The interaction 
drastically changes the evolution of the cc distribution 
in position space, and this will be demonstrated again in 
the full numerical simulation of the next section. 
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FIG. 4: (Color online.) Distribution over quark pair separa- 
tion at fixed T — l.STc after 9 fm/c, with (red squares) and 
without (green triangles) the cc potential. 



FIG. 5: Energy distribution of the cc pairs (in the center of 
mass frame of the pair) after Langevin evolution at a fixed 
temperature T = 1.05Tc for a time t — 9 fm/c long enough 
for quasiequilibrium to be reached. 



III. LANGEVIN EVOLUTION OF THE cc 
PAIRS 

A. Langevin evolution in static medium: 
quasiequilibrium 

Before we turn to expanding fireball, wc first study the 
evolution of cc pairs at fixed T = 1.5Tc and in the absence 
of hydrodynamical flows. The first thing we would like to 
demonstrate is the strong influence of the heavy quark in- 
teraction. The resulting distribution over interquark sep- 
aration at time 10 fm/c are shown in Fig. IIII Al with the 
interaction (red squares) and without (blue triangles). 
The value for a given pair's separation r is weighted by 
^ so that the spatial phase space of the distribution is 
divided out. With the interaction on, we find the same 
behavior of a "slowly dissolving lump" as that seen in 
the solution of the Fokker-Planck equation when the in- 
teraction is "turned on" . 

Further study of this has shown convergence of its 
shape to a particular one, which persist for a long time 
and which we would call "quasiequilibrium"'^. While the 
true equilibrium of course corresponds to complete dis- 
solution of a single cc pair, it turns out that leakage to 
large distances affects the distributions of separation and 
energy in normalization only. In Fig. IIII Al is the energy 
distribution for the ensemble of pairs at r = 9 fm/c after 
evolving under Langevin dynamics at a fixed temperature 
1.05Tc, and it is shown to be the same distribution, up to 
normalization and statistical uncertainty, as the distribu- 
tion reached by the pairs in a full heavy-ion simulation 
of the most central collisions. 



^ This situation should not be confused with stationary nonequi- 
librium solutions of the Foklcer-Planck equation, in wliich there 
is a constant flow through the distribution because of matching 
source and sink. 



We show the energy distribution in this region because 
it is related to a very important issue of the charmonium 
production, namely production of -0 , x states and subse- 
quent feeddown into the J/ip- When a quasiequilibrium 
distribution is reached, the production ratios of charmo- 
nium states are stabilized at thermal (statistical) model 
values, in spite of the fact that the overall normalization 
continues to shrink due to leakage into infinitely large 
phase space at large distances. 

(The energy distribution itself contains a Boltzmann 
factor but also the density of states. A model case of 
purely Coulomb interaction allows one to calculate it in 
a simple way: as shown in Appendix |^ we found that 
in this case the absolute shape of the quasiequilibrium 
distribution is reproduced as well.) 

The existence of quasiequilibrium is in good correspon- 
dence to observations. It was noticed already a decade 
ago [i^l that the SPS data on centrality dependence of 
N^i/Nj/^ ratio approached the equilibrium value (for 
chemical freezout) 

^ = exp{~/^MlT,n) (24) 

with the chemical freezout at Tch = 170 MeV, as is ob- 
served for ratios of other hadronic species. 

One possible explanation of it can be charmonium 
recombination (from independent charm quarks) at 
chemical freezout, advocated by [3| and others. How- 
ever our findings show that the same ratio naturally ap- 
pears even for a single cc pair dissolving in a thermal 
medium, in a "quasiequilibrium" occurring at the leakage 
stage. Especially at SPS, when statistical recombination 
requires a charm density which is too large, this is an 
attractive possibility. 
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FIG. 6: (Color online.) Evolving energy distribution for an 
ensemble of cc pairs at time moments t = 2, 3, 10 fm/c (cir- 
cles, squares and triangles, respectively). 



C. Langevin motion of cc pair in an expanding 
fireball 

Finally, we model the motion of a charm quark pair in 
an evolving fireball. We use the same framework and pro- 
grams used in [2^ to examine motion of a single charm 
quark for propagation of an interacting pair. 

We start with large number of cc produced with 
PYTHIA pQCD event generation, and randomly place 
them in position space throughout the fireball, using a 
Monte-Carlo Glauber calculation. 

Then, the pairs are evolved in time according to the 
Langevin equations: 

5 = + e - VC/ (26) 



B. Production of the initial cc pairs 

We start with cc events produced with PYTHIA, a 
particle physics event generator [s^. PYTHIA yields 
cc pairs through a set of perturbative QCD events: 
through Monte-Carlo it will select initial partons from 
the parton distribution function of one's choosing and 
proceed through leading-order QCD differential cross sec- 
tions to a random final state. The pt and rapidity distri- 
butions of charm produced in pp collisions is believed to 
be rather adequately represented. 

By using PYTHIA we do not however imply that it 
solved many open issues related with charm production, 
such as color octet versus singlet state. It also leaves open 
the very important issue of feeddown from charmonium 
excited states (see below). One more open question - 
needed to start our simulations ~ is how to sample the 
distribution in position space. Indeed, each pQCD event 
generated in PYTHIA is a momentum eigenstate without 
any width, so by Heisenberg's uncertainty relation they 
are spatially delocalized, which is unrealistic. 

We proceed assuming the form for the initial phase- 
space distribution to be 

Pinit(r, p) « Ppyth(p) exp(-rV2a2) (25) 

By setting a = 0.3 fm one can tune energy distribution 
to give a reasonable probability for the formation of Jf-i/j 
state in pp events. 

However this does not yield correct formation probabil- 
ities of X, ■ It is hardly surprising, since for example the 
Tp wavefunction has a sign change at certain distances, 
and the exact production amplitude is required for a pro- 
jection, an order-of-magnitude size estimate is not good 
enough. Since feeddown from those states contributes 
about 40% of the observed J/V', we simply refrain from 
any quantitative statements about pp (and very periph- 
eral AA) collisions, focusing only on distributions after 
some amount of time spent in QGP. 



where ^ corresponds to a random force and rj the drag 
coefficient. The condition that the Langevin equations 
evolve a distribution towards thermal equilibrium gives 
the relation = 2MTr]5ij5{t - t'). We pro- 

ceed here with our diffusion constant set by the results 

of [2^ to be rjD = yEm- ^'^'^ discussed earlier with 
our potential as V instead of F. 

Now we examine the evolution of the quark pairs as 
discussed before, examining pairs at mid-rapidity in a 
boost-invariant, 2-dimensional ultra-relativistic gas simu- 
lation, the same hydrodynamical simulation used in [2^ . 
We stop the Langevin-with-interaction evolution when 
T < Tc- The distribution over energy at different mo- 
ments of the time is shown in Fig llll Bl 

We will discuss our results subsequently, at two dif- 
ferent levels of sophistication. First we will show results 
with only the total number of bound states monitored, 
and then we will show results where the different char- 
monia states and feeddown to J/iJj are considered and 
compare these results with PHENIX data. 

How we determine whether or not to call a cc pair 
in our simulation bound, and what particular charmo- 
nium state the pair exists as if it is bound, is discussed 
in Appendix IB] Fig. IIII CI shows the number of "bound" 
cc pairs as a fraction of their total number, monitored 
throughout the course of the simulation as a function of 
time. Note that realistic hydrodynamical/Langevin sim- 
ulation agrees qualitatively with the analytic results of 
Section ini During the first fm/c one finds some boosts 
in the probability for a cc pair to be bound, due to rapid 
thermalization in momentum space. Later the probablity 
falls due to the slow diffusion in position space. This fig- 
ure emphasizes our main qualitative finding, the survival 
probability of J/ tp on the order of a half. 
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FIG. 7; (Color online.) Probability of cc pairs to be bound 
at RHIC Au+Au, y/s = 200 GeV, mid-rapidity. 



D. Shadowing and "normal" absorption 

Experimental data include not only the "QGP sup- 
pression" we study in this work, but also (i) the initial- 
state effects (modified parton distributions in "cold nu- 
clear matter" ) plus (ii) the so called "normal nuclear ab- 
sorption" . The way we have chosen to display PHENIX 
data [4l| is as follows: before we compare those with our 
results we "factor out" the cold nuclear matter effects, 
by defining (for any given rapidity y) the following ratio 
of Au-|-Au and d-|-Au data 



R 



anomalous 



AA 



(y) 



nPHENIX 



(y) 



RpA{y)RpA{-y) 



(28) 



to be called the "anomalous suppression" . In principle 
those include only data: but unfortunately the large dAu 
sample taken in 2008 is not yet analyzed (at the time of 
this writing) , while the 2003 set has error bars which are 
too large. This forces us to use a model at this point, fol- 
lowing Kharzeev et aljsl] with RpA = exp{—aabs{L)no), 
where (i) is the mean path length of the J/ ip through nu- 
clear matter, hq is the nuclear density, and Cabs is the nu- 
clear absorption cross-section (parametrized from 39] to 
be O.lfm^ for rapidity y = 0). Finally, for rapidity y — 0, 
we rewrite this as {RpA{v = 0))^ = exp(-cTahs (Npart), 
where (Npart) is the density per unit area of participants 
in the collision plane. We further used Glauber model 
calculations 40] to determine (Npart) for a given A^part 
at PHENIX. We divide each Au+Au data point from 
PHENIX by this quantity and caU it R'Xa"^"^"'"' (v = 0) 
plotted as points in Figs lIIIDl and (TUl to be compared 
with our simulation. 



E. tp production and feeddown 



'/Nj 



in our simula- 



Next we calculate the ratio N, 
tions, for different centralities 
NA38/50/60 measurements of this ratio at the SPS, but 



There are well known 



FIG. 8: (color online) Points show the magnitude of the 
anomalous suppression at mid-rapidity _RJ^'^'""'°"*(j/ = 0) 
versus the centrality (the number of participants), using 
Au+Au PHENIX data. The curve is the probability to be 
bound (determined by the energy projection) at the end of 
the QGP era, when T = T^. 



at RHIC it has been measured so far only in pp col- 
lisions by the PHENIX detector [H] to be 0.14 ± 0.04, 
which makes the ratio of direct ip to J/ip particles 0.24 as 
in 43]. Hopefully higher luminosity at RHIC will make 
possible a future measurement of this ratio in Au+Au 
collisions of various centralities. 

We calculate A^^' /Nj/^ as follows: (i) first we run our 
simulation and determine the distributions f{E) over the 
cc pair energy E = Ecu — '2.Mc (in the pair center of 
mass frame); (ii) then we compare those to quasiequilib- 
rium ones from simulations at fixed temperature (slightly 
above Tc) fo{E). Since both are done for the same in- 
teraction, in the ratio f{E)/fo{E) the density of states 
drops out. This ratio tells us how different the actual 
distribution is from that in quasiequilibrium. (iii) Then 
we form the double ratios, at two relative energies corre- 
sponding to ip and J/ip masses (minus 2mcharm) 



R 



/(.8 GeV) , /(.3 GeV) 



/- 



/o(.8 GeV)' /o(.3 GeV) 



(29) 



This now includes nonequilibrium effects for both of 
them. Finally (iv) we switch from continuum clas- 
sical distributions to quantum one, assuming that in 
quasiequilbrium the relation ((24)) holds. If so, the particle 
ratio is a combination of nonequilibrium and equilibrium 
factors 



N„ 



N 



^^'/V.exp(-Am/T) 



(30) 



The double ratio ( or exp{AM /T)N^, /Nj/^) is plotted 
vs centrality in Fig. [HI As one can see, it goes to unity 
for the most central collisions: so quasiequilibrium is 
actually reached in this case. For mid-central bins the ip 
production is about twice larger because of insufficient 
time. This is to be compared to the experimental pp 
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FIG. 9: (Color online.) The double ratio R^' defined in 
(|29p versus centrality (number of participants). One point 
(green box) at Npart = 2 corresponds to experimental data 
for and the direct J/ip , for pp collisions. 
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FIG. 10: (Color online 

-jartomalous 



are our model, with (solid, 
lower) feeddown. 



upper) and without (dashed, 



value for the ratio, which is about 5. (We remind 
the reader that PYTHIA plus our classical projection 
method does not work for pp collisions.) 

Finally, we use this result to estimate the effect of feed- 
down from higher states. To do this, we write the final 
number of J/ip particles observed as the number of di- 
rectly produced J/tp particles plus the number of J/tp 
particles produced from feeddown from higher charmo- 
nium states: 



N 



final 



-\jdirect 



1 + 



5i ^M, 
y)exp( 



31) 



where i is summed over the xi, X2, and ipi particles which 
contribute significantly to feeddown, Bi represents their 
branching ratio into J/ip + and gi is the degeneracy 
of the state (divided by 3, the degeneracy of the J/i/')j 
Ami is the mass difference between the i-th state and the 
J/^. The R^' 1^ is the non-equilibrium factor discussed 
above: it is factored outside the sum because it is very 
similar for all these states. 

Now we are ready to discuss centrality dependence of 
the J jip production including the feeddown. We define 
for each centrality direct N'j;J'^''^{b) as the total number 
of cc pairs in our ensemble with energy (in its rest frame) 
less than E < 2Mcharm + 0.33 GeV. The feedown gets 
its dependence on centrality from R^,, /,^{b) determined 
from simulation. 

The absolute normalization of the results deserves spe- 
cial discussion. We predict the absolute probability of 
J/V' production, both direct and with feeddown, nor- 
malized per cc pair produced in the same collisions (e.g. 
centrality bin). Unfortunately, the total cross section of 
charm production is not yet measured with sufficient ac- 
curacy to normalize results this way. 

The usual way to present these results is in the form of 
the so called Raa ratio, relating production in Au+Au 



at given centrality to that in pp, times the theoretical 
(Glauber) predictions for the number of hard collisions. 
In other words, the parton model for cc is used. Unfor- 
tunately, the experimental data about feeddown in pp 
are still uncertain enough to produce sufficiently large 
scale ambiguity. We cannot obtain this from our the- 
ory as well because (as explained above) our classical 
projections do not work adequately for pp (and very pe- 
ripheral) collisions. Thus we see at the moment no suf- 
ficiently accurate way of absolute comparison with the 
data. Because of that, we simply normalize our results 
assuming that there is no QGP suppression ( meaning 
j^anomaious^j^^^^^ < 100) = 1). The rcsults normalized 
like this are shown in Fig. [TOl we conclude that although 
feeddown is not large, taking it into account helps bring 
the shape of the anomalous suppression closer to obser- 
vations. 



F. Including the "mixed phase" 

In our work so far, we have only examined the evo- 
lution of the cc pairs during the QGP phase, stopping 
the evolution wherever the fluid's temperature reached 
Tc- However, in order to understand the evolution of 
charmonia to their hadronization we need to model the 
dynamics of the charm quarks also through the "mixed 
phase" , also known as the near-Tc region. 

In various hydrodynamic models which describe heavy 
ion collisions, the region of roughly T = 0.9 — 1.1 Tc is 
treated as a separate "mixed phase" distinct from QGP 
and hadronic phases. Indeed, it has a very different equa- 
tion of state: while the temperature and pressure remain 
nearly constant, the energy and entropy densities jump 



11 



by a large factor ^. 

What is very important for the present paper is that 
the near-Tc region occupies a significant fraction of the 
fireball evolution, in spite of being a very narrow interval 
in terms of T . Indeed, one should rather focus not on 
T but on the entropy density s, which shows a simple 
monotonous decrease with time s ^ 1/r for all three 
phases. 

For a quantitative description of the mixed phase we 
used hydrodynamical calculations, known to describe ra- 
dial and elliptic flow well, such as the work by Shuryak 
and Teaney [i^. It follows from their solution that the 
"mixed phase" persists for about 5 fm/c after the decon- 
fined phase, which is comparable to the lifetime of the 
deconfined phase at the very center of the fireball. Thus 
it is by no means a small effect, and should be included 
in any realistic treatment of a heavy-ion collision. 

The flow during this time was found to be well ap- 
proximated by a Hubble-like expansion with radial ve- 
locity V = Hr and time-independent H w 0.075 fm^^ 
for central collisions. For a collision with a nonzero im- 
pact parameter (below we will consider h = 6.5 fm), the 
anisotropy of this expansion can be parameterized simi- 
larly: 

Vi = HiXi (32) 

with i ^ 1,2 and i?^ = 0.078 fm"\i/j, = 0.058 fm"^: 
thus anisotropy is only about 80% by this late stage. It 
is fair to say that we have a fairly reasonable understand- 
ing of how the medium flows for these later stages: thus 
in our simulations we have used those parameterizations 
instead of numerical solutions to hydrodynamics, which 
were necessary for the QGP phase. 

Let us start with two extreme scenarios for the dynam- 
ics of the charm quarks during this phase of the collision: 

1. ) the charm quarks are completely "stopped" in the 
medium, so that they experience the same Hubble-like 
flow as matter; 

2. ) cc pairs do not interact at all with the medium near- 
Tc, moving ballistically with constant velocity for the cor- 
responding time in the collision. 

If the first scenario were true, the effect of Hubble 
flow would be to increase all momenta of particles by the 
same multiplicative factors Pi{t) = Pi(0) exp{Hit). With 
sufficiently high drag, Langevin dynamics would bring 



Although the exact nature of matter in the near-Tc region is not 
yet understood, let us mention that the original "mixed phase" 
description, originating from the notion of the first-order phase 
transition, cannot be true, as "supercooling" and bubble for- 
mation expected in this case are not observed experimentally. 
Lattice gauge theory suggests a smooth crossover-type transi- 
tion, with a high but finite peak in the specific heat. Recently 
there has been renewed interest in this region, after the so-called 
"magnetic scenario" for it has been proposed |46l |47|| , describ- 
ing it as a plasma containing a significant fraction of magnetic 
monopoles. 
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FIG. 11: (Color online.) The charm pt-distribution after the 
mixed phase compared with the distribution without flow, the 
distribution orignating from PYTHIA, and the distribution 
before the mixed phase 



the charm quarks rapidly to a thermal distribution, and 
since M > > T it is a good approximation in this case to 
say that the heavy quarks have been " stopped" . How- 
ever, we will show below that at the "realistic" value 
used for the drag ijc this does not happen during the 
time allocated to the mixed phase, there is instead on- 
going "stopping" of the charm quarks relative to fluid 
elements. (This also will be important for the evolution 
of the azimuthal anisotropy V2{pt) for single charm and 
for charmonium). The second scenario predicts V2{pt) 
for single charm quarks which is far smaller than what 
is measured. We do not consider this scenario further 
even though something might be said for modelling char- 
monium in the mixed phase as interacting far less than 
single charm. 

Several single charm pt-distributions are shown in Fig- 
ure [11] (normalized for simplicity to unity). The initial 
distribution after hard production predicted by PYTHIA 
is the largest at high pt: this is compared to the Langevin 
evolution before (squares) and after (triangles) the mixed 
phase, for a semi-central RHIC collision (6 = 7 fm). 
In order to see that radial flow still is important, we 
have also shown what happens if Langevin evolution hap- 
pen on top of unmoving fixed-T plasma (circles). This 
comparison demonstrates once again the main point of 
this paper, that for charm quarks and charmonium in a 
heavy-ion collision equilibration is never complete, even 
in momentum space: so the specific timescales of differ- 
ent phases of matter are of fundamental importance. 

Unfortunately, in the near-Tc region it is much less 
clear how to describe the c — c interaction. As we have 
learned from lattice data, the difference between free- 
energy and potential-energy potentials are very drastic 
in this case: in the former case the tension (the force 
from the linear potential) disappears while in the latter 
it becomes about 5 times stronger than it is in vacuum. 
As discussed in refs[4^|4^, the latter is presumably due 
to a metastable electric flux tube. 
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FIG. 12: (Color online.) The azimuthal anisotropy versus 
transverse momentum for both single charm and for J/tp. 



Which potential to use depends on timescales of the 
c — c dynamics, which are not well understood at this 
point. Therefore we took for now a conservative ap- 
proach, assuming that at the near-Tc stage charm pairs 
interact according to the simple Coulomb interaction 
V = —OLsjr. Additionally, in our model for this phase 
we assume the interaction of the charm quarks with the 
medium can be modelled with the same Langevin dy- 
namics with the temperature approximated as a fixed 
T = Tc and the flow given as above. We found that with 
the simple Coulomb potential used in the mixed phase, 
the survival probability dropped slightly but not signifi- 
cantly: and although we do not discuss other possibilities 
in this work further, in principle this can be changed if 
the potential to be used has significant tension. 

One final interesting observable would be a measure- 
ment of charmonium elliptic flow, characterized by the 
azimuthal anisotropy parameter V2 —< cos{2<j)) >, in- 
duced by ellipticity of the flow on charmed quarks. A 
measurement with low statistics has been already made 
at PHENIX [13]: both PHENIX and STAR are working 
on higher statistics data on tape now. The result of our 
calculation of V2{pt)^ both for single charm quarks and 
for the J/'0, is shown in Figure fT2l 

Greco, Ko, and Rapp also made predictions for the 
V2 for D and J/ip i51|, based on a completely different 
scenario: in their case charm distributions are completely 
equilibrated and the charmonium states coalesce from 
them at hadronization. In spite of such difference, our 
predictions are similar: V2 of J /ijj should be less than the 
V2 of single charm for low pt, but then increase past the 
V2 of single charm a.i pt > i GeV . This shows, that the 
observation of charmonium V2 can not be considered the 
argument for coalescence. 



the sQGP stage using hydrodynamics-|-Langevin. The 
main elements of the paper are: (i) inclusion of the inter- 
action force between charm quarks and (ii) emphasis on 
deviations from equilibrium during the finite QGP life- 
time. 

Our main finding is that the lifetime of sQGP is not 
sufficient to reach the equilibrium distribution of the 
pairs in space, allowing for a significant fraction of J /ill 
^ 1/2 to survive. This probability for charmonium dis- 
sociation in sQGP is larger than in earlier perturbative 
estimates, or for Langevin diffusion ignoring mutual in- 
teraction. 

That is why there is no large difference between sup- 
pression at RHIC and SPS, in spite of the longer QGP 
lifetime in the former case. While the momentum re- 
laxation is rather rapid, we found that later evolution 
reaches the so-called quasiequilibrium regime, which is 
maintained during all time of QGP expansion. The spa- 
tial distributions after some time develop a "core" in 
which cc pairs remain in close proximity due to the re- 
maining effective attraction in sQGP, combined with a 
relatively slow leakage into a spreading tail toward large 
r. We propose quasiequilibrium as the clue to an expla- 
nation of the apparently thermal ratios of ■0 and J /'lp^ 
especially at SPS. 

The shape of the centrality dependence of our survival 
probabality is in agreement with data. Therefore, al- 
though we have not yet directly evaluated "nondiagonal" 
recombination, we think that most of the J /^p observed at 
SPS and RHIC are still from the "diagonal" pairs. This 
and other issues will of course be clarified, as more simu- 
lations and data in different kinematic domains become 
available. 

The calculation is extended to the near-Tc region - 
known as a "mixed phase" in hydrodynamical calcula- 
tions. Its duration for RHIC collisions is about 5 fm/c, 
comparable to that of the QGP. We have used simple 
Hubble-like parameterization and minimal Coulomb po- 
tential, and predict both the charmonium pt spectra and 
azimuthal asymmetry parameter V2 which is expected to 
be measured soon . 
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IV. SUMMARY 



APPENDIX A: CLASSICAL VS. QUANTUM 
DYNAMICS 



We have studied a relaxation process of a cc pair pro- 
duced by hard processes in heavy ion collision throughout 



In this paper, we take phase-space distributions of 
cc pairs and evolve them according to the Fokker- 



13 



Planck/Langevin equations, which describe nonequihb- 
rium evolution of phase-space distributions during the 
QGP era. After it is finished and the medium returns 
to the hadronic phase, our classical distribution f{x,p, t) 
has to be projected into charmonium quantum states. 

He we examine how classical and quantum dynamics 
correspond to each other in equilibrium, when both are 
easily available. We simplify by examining the thermal 
distributions for a Coulombic system, with V ^ 1/r. One 
can calculate the density of states for the classical system: 

Z = [drdp{4TTfr'^p'^exp{-{p^/2fi-e^/r)/T) 



dE I drdp{ATTfr^p^exp{-E/T)S{E-p^/2^i+l/r) 



x'Indt 0.9188/4 
Constant 1 .777e-07 ± 8.622e-08 
Power -2.087 ±0.1 33 



10-*- 




-0.1 
Binding Energy [GeV] 



= / dEcxp{-E/T) [ 

J Jrn 



dp 



2 9 



ipy2^i-Er 



/ m,ax{0,\/2jrE 

As one can see, this integral is divergent for E > 0. This 
means that this distribution is not normalizable and in 
thermal equilibrium all pairs are ionized. For i? < 0, we 
see from examining the partition function 

p{E) (X (-^)-5/2 (Al) 

Now we calculate the quantum mechanical density of 
states for _E < 0, 



piE) = ^^^,n'S{E+^) 



(A2) 



which can be approximated by considering an integral: 

E ^ 

dE'p{E') = j:°g^t^9{E + -) 



~ E "^1"^ for large enough E 
Thus, p(E) - £^-5/2 fQj. ^ ^jQgg 

zero: therefore for 

highly excited states the correspondence principle holds 
and the classical thermal distribution is recovered. How- 
ever, classical density of states is not so good of an ap- 
proximation to the quantum- mechanical density of states 
for the lowest bound state. 

We also tested whether Langevin simulations leads to 
the correct classical density of states, after some relax- 
ation time. The result of evolving an ensemble of cc pairs 
at a fixed temperature according to classical Langevin dy- 
namics shown in Fig |13l shows that equilibrium is indeed 
obtained. 



APPENDIX B: COALESCENCE PROBABILITY 

After all distributions at the end of the QGP era are 
determined, the next step is to calculate the probabilities 
of it materializing as a particular charmonium state. 

One approach is based on the Wigncr probability dis- 
tribution. For any wavefunction '(/'(r), 

W^V-(r,p) = ^ /d'#*(r + |Mr-|)e'P-^ (Bl) 



FIG. 13: The points show the density of states for an ensemble 
of cc pairs at a fixed temperature T = l.STc obtained by long- 
time Langevin evolution, compared to a fit with power 2.1 
(curve), close to classical nonrelativistic thermal distribution. 



The wave function for J/i/; was easily calculated numer- 
ically and then fitted to a Gaussian '^i{x) oc e~^^^'^ . 
This leads immediately to the following parametrization 
of the J/i/i's Wigner distribution's (here the constant is 
in fm^^) : 

Wji^ ix,p)^ ^exp{-7.6r^ - p^l.e) (B2) 

Many properties for the Wigner transform so defined 
may be discussed but for our purposes here we should 
note that 



W^{x,p)dp=\^{x)\' 



(B3) 



and for another wavefunction xi^y^ Wigncr transform 
W^ix,p) 



■^ix,p)W^{x,p)dxdy (B4) 



because this is what we need to properly normalize our 
distributions and to compute overlaps. 

So, for our phase-space distribution at any given time 
/(x, p, t), we model our probability of pair being mea- 
sured in the J/V'-state as 

Pj/^ (r) = (2^)3 J /(r, p, t)Wj/^ (r, p)^^^^^^ (B5) 

We project the pairs onto the J/i/i-state not only at the 
onset of the calculation, but also throughout the time 
evolution, monitoring in this way the probability of J/tp 
production. We use this approach to estimate the coa- 
lescence probability for the distribution in Section iHl 

Later on we switched to the "energy projection 
method" , which is ultimately used for projection into 
X, ip states as well as into J/ip . It is based on the distri- 
bution over cc energy, in the pair rest frame, calculated 
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with the (zero temperature) Cornell potential. One pro- 
jection, used in Fig llll Dl is to all bound states, defined 

by 

EcM = Vcorneii i^) + Jf < ^'^^ GeV. (B6) 
Later on, we differentiate between various charmonium 



states again by examining the cc pair's energy in its rest 
frame, and comparing it with the energies of solutions to 
Schrodinger's equation using the Cornell potential. For 
example, the lowest s-state solution, using the charm 
quark mass, has energy E = 0.33 GeV, therefore we 
count all cc pairs in our simulation with energy below 
0.33 GeV as existing in the J/tp state. 
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